

icr_rands::icr_rands(void){
  //constructor
  initialize_random_numbers();
}
//////////////////////////////////////

void icr_rands::initialize_random_numbers(void){
  // ---
  //a bunch of stuff has to be set up for the seed to work.
  // First we get the time
  time_t rawtime;
  struct tm * timeinfo;
  time (&rawtime);
  timeinfo=localtime(&rawtime);
  cout<<"ICR: Start at time: "<<asctime(timeinfo)<<endl;
  //now we define the seed with microsecond precision
  timeval rawtime2;
  gettimeofday(&rawtime2, NULL);
  //This has defined the time on the computer which we use as a seed            
  long t1=rawtime2.tv_sec + rawtime2.tv_usec;
  if(seed != 0) t1=seed;

  //using seed is the pointer to the actual seed
  using_seed = new base_generator_type (t1);
  cout<<"ICR: Seed used "<<t1<<endl;

  boost::uniform_real<> uni_dist(0,1);
  
  uni=new  boost::variate_generator<base_generator_type, boost::uniform_real<> >                (*using_seed, uni_dist);

for(int i=0;i<10000;i++) (*uni)(); //add some entropy to seed

 boost::normal_distribution<> pdf(0.0,1.0);
 normgen = new normal_generator (*using_seed, pdf);


}
//////////////////////////////////////
gamma_generator icr_rands::rgamma_gen(double alpha){
  boost::gamma_distribution<> pdf(alpha);
  gamma_generator generator(*using_seed, pdf);
  return generator;
}
//////////////////////////////////////

rmvnorm icr_rands::mvnorm_gen(rowvec mu,mat sigma){
   return rmvnorm(mu, sigma, this);
}

//////////////////////////////////////
riwish icr_rands::riwish_gen(double nu,mat S){
   return riwish(nu, S, this);
}
//////////////////////////////////////
rowvec icr_rands::rdirichlet(rowvec alphas){
    //draw from N gamma deviates each with the proper alpha, then divide by sum
    double sum=0;
    double drawn;
    rowvec output(alphas.n_elem);
    for(int i=0;i<(int) alphas.n_elem;i++){
       drawn  = rgamma_gen(alphas[i])();
       sum += drawn;
       output(i) = drawn;
    }
    for(int i=0;i<(int) alphas.n_elem;i++){
       output(i) /= sum;
    }
    return output;
}

//////////////////////////////////////
long icr_rands::random_index(rowvec cumsum0){
   //rowvec cumsum0 = cumsum(wgts);
   double rand0 = rand01();
   for(int i = 0;i<cumsum0.n_elem;i++){
      if(rand0 <= cumsum0[i]) return i;
   }
   cumsum0.print("cumsum");
   cout<<rand0<<endl;
   cout<<"icr_rands::random_index Shouldn't get here!!!"<<endl;
   exit(1);
   return cumsum0.n_elem-1;
}


//////////////////////////////////////

icr_rands::~icr_rands(void){
  //destructor
  delete uni;
}
